# Cleaning the Work Environment ####
rm(list=ls(all=T))
gc()

getwd()

library(sjPlot)
library(multcomp)
library(stargazer)

#### Loading the Dataset ####

s2 <- readRDS(file = "SEEBSS - Study 2.rds")

#### Data Preperation ####

table(s2$Q31_1, exclude = NULL)
s2$Q31_1[s2$Q31_1==99] <- NA
names(s2)[names(s2)=="Q31_1"] <- "Self_Conservatism"
table(s2$Self_Conservatism, exclude = NULL)

table(s2$Q36, exclude = NULL)
s2$Q36[s2$Q36==99] <- NA
s2$Q36[s2$Q36==2] <- 0
names(s2)[names(s2)=="Q36"] <- "sex" # 1 - female, 0 - male
table(s2$sex, exclude = NULL)
s2$sex <- 1 + s2$sex*-1 # recoding # 0 - female, 1 - male
table(s2$sex, exclude = NULL) 

s2$Q38[s2$Q38==99] <- NA
names(s2)[names(s2)=="Q38"] <- "education"
table(s2$education, exclude = NULL)

s2$Q7_1[s2$Q7_1==99] <- NA
names(s2)[names(s2)=="Q7_1"] <- "self_left_right"
table(s2$self_left_right, exclude = NULL)

s2$Affiliation <- NA
s2$Affiliation[s2$Q3==99] <- 0
s2$Affiliation[s2$Q3==1] <- 2
s2$Affiliation[s2$Q3==2] <- 1
table(s2$Affiliation, exclude = NULL)

table(s2$Q4_TEXT, exclude = NULL)

table(s2$Q4, exclude = NULL)
s2$Partisanship <- NA
s2$Partisanship[s2$Q4==1] <- "Justice and Development Party (AKP)"
s2$Partisanship[s2$Q4==2] <- "Republican People's Party (CHP)"
s2$Partisanship[s2$Q4==3] <- "Other Parties on the Left"
s2$Partisanship[s2$Q4==4] <- "Other Parties on the Right"
s2$Partisanship[s2$Q4==5] <- "Other Parties on the Right"
s2$Partisanship[s2$Q4==6] <- "Other Parties on the Right"
s2$Partisanship[s2$Q4==7] <- "Other Parties on the Left"
s2$Partisanship[s2$Q4_TEXT == "Emek Partisi" | 
                  s2$Q4_TEXT == "EMEP" | 
                  s2$Q4_TEXT == "ÌÐDP" | 
                  s2$Q4_TEXT == "TÌ_rkiye €¡\u0081ÙÌ¤i Partisi" | 
                  s2$Q4_TEXT == "TÌ_rkiye komÌ_nist partisi" |
                  s2$Q4_TEXT == "TÌ_rkiye KomÌ_nist Partisi" |
                  s2$Q4_TEXT =="TÌ_rkiye Kominist partisi" |
                  s2$Q4_TEXT =="TÌ_rkiye komunist partisi" |
                  s2$Q4_TEXT =="Tip" |
                  s2$Q4_TEXT =="Tkp" |
                  s2$Q4_TEXT =="TKP"] <- "Other Parties on the Left"
s2$Partisanship[s2$Q4_TEXT == "Cumhur ittifak€±" | 
                  s2$Q4_TEXT == "Cumhur Ittifak€±" |
                  s2$Q4_TEXT == "ÌÐtÌ_ken Birli€Ùi Partisi" |
                  s2$Q4_TEXT == "Liberal Demokrat Parti" |
                  s2$Q4_TEXT =="Sa€Ùduyu"] <- "Other Parties on the Right"
s2$Partisanship <- ifelse(is.na(s2$Partisanship),"Bipartisan or Non-Responsive", s2$Partisanship)
table(s2$Partisanship, exclude=NULL)
s2$Partisanship <- factor(s2$Partisanship, ordered = F) 
table(s2$Partisanship, exclude=NULL)

s2$bipartisan_na <- NA
s2$bipartisan_na[s2$Partisanship=="Bipartisan or Non-Responsive"] <- 1
s2$bipartisan_na[s2$Partisanship!="Bipartisan or Non-Responsive"] <- 0
table(s2$bipartisan_na, exclude = NULL)

s2$akp <- NA
s2$akp[s2$Partisanship=="Justice and Development Party (AKP)"] <- 1
s2$akp[s2$Partisanship!="Justice and Development Party (AKP)"] <- 0
table(s2$akp, exclude = NULL)

s2$chp <- NA
s2$chp[s2$Partisanship=="Republican People's Party (CHP)"] <- 1
s2$chp[s2$Partisanship!="Republican People's Party (CHP)"] <- 0
table(s2$chp, exclude = NULL)

s2$other_right <- NA
s2$other_right[s2$Partisanship=="Other Parties on the Right"] <- 1
s2$other_right[s2$Partisanship!="Other Parties on the Right"] <- 0
table(s2$other_right, exclude = NULL)

s2$other_left <- NA
s2$other_left[s2$Partisanship=="Other Parties on the Left"] <- 1
s2$other_left[s2$Partisanship!="Other Parties on the Left"] <- 0
table(s2$other_left, exclude = NULL)

mean(s2$akp, na.rm = T)
s2$beta_akp <- (s2$akp - mean(s2$akp, na.rm = T))
table(s2$beta_akp, exclude = NULL)

mean(s2$chp, na.rm = T)
s2$beta_chp <- (s2$chp - mean(s2$chp, na.rm = T))
table(s2$beta_chp, exclude = NULL)

mean(s2$other_right, na.rm = T)
s2$beta_other_right <- (s2$other_right - mean(s2$other_right, na.rm = T))
table(s2$beta_other_right, exclude = NULL)

mean(s2$other_left, na.rm = T)
s2$beta_other_left <- (s2$other_left - mean(s2$other_left, na.rm = T))
table(s2$beta_other_left, exclude = NULL)

table(s2$Group, exclude = NULL)

s2$Cancer_Prime <- NA
s2$Cancer_Prime[s2$Group == "Treatment 1"] <- 1
s2$Cancer_Prime <- ifelse(is.na(s2$Cancer_Prime), 0, s2$Cancer_Prime)
table(s2$Cancer_Prime, exclude = NULL)

s2$Terror_Prime <- NA
s2$Terror_Prime[s2$Group == "Treatment 2"] <- 1
s2$Terror_Prime <- ifelse(is.na(s2$Terror_Prime), 0, s2$Terror_Prime)
table(s2$Terror_Prime, exclude = NULL)

s2$General_Mortality_Prime <- NA
s2$General_Mortality_Prime[s2$Group == "Treatment 3"] <- 1
s2$General_Mortality_Prime <- ifelse(is.na(s2$General_Mortality_Prime), 0, s2$General_Mortality_Prime)
table(s2$General_Mortality_Prime, exclude = NULL)

table(s2$Q14_7, exclude = NULL)
s2$Q14_7 <- ifelse(s2$Q14_7>11, NA, s2$Q14_7)
s2$Q14_7[s2$Group=="Control"] <- 0

s2$cancer_prime_fear <- s2$Q14_7
table(s2$cancer_prime_fear, exclude = NULL)

table(s2$Q15_7, exclude = NULL)
s2$Q15_7 <- ifelse(s2$Q15_7>11, NA, s2$Q15_7)
s2$Q15_7[s2$Group=="Control"] <- 0

s2$terror_prime_fear <- s2$Q15_7
table(s2$terror_prime_fear, exclude = NULL)

table(s2$Q16_7, exclude = NULL)
s2$Q16_7 <- ifelse(s2$Q16_7>11, NA, s2$Q16_7)
s2$Q16_7[s2$Group=="Control"] <- 0

s2$general_mortality_prime_fear <- s2$Q16_7
table(s2$general_mortality_prime_fear, exclude = NULL)

s2$priming_fear <- rowSums(s2[, c("cancer_prime_fear", "terror_prime_fear", "general_mortality_prime_fear")], na.rm = T)
table(s2$priming_fear, exclude = NULL)
s2$priming_fear <- ifelse(is.na(s2$cancer_prime_fear) & is.na(s2$terror_prime_fear) & is.na(s2$general_mortality_prime_fear), NA, s2$priming_fear)

#### Manipulation Check ####

s2$Q17_1[s2$Q17_1==99] <- NA
s2$Q17_2[s2$Q17_2==99] <- NA

s2$MC1 <- NA
s2$MC1[s2$Q17_1==6 & s2$Q17_2==8] <- 1
s2$MC1[(s2$Q17_1!=6 | s2$Q17_2!=8) & (!is.na(s2$Q17_1) & !is.na(s2$Q17_2))] <- 0
s2$MC1 <- ifelse(is.na(s2$MC1), 0, s2$MC1)
table(s2$MC1, exclude = NULL)

s2$Q18_1[s2$Q18_1==99] <- NA
s2$Q18_2[s2$Q18_2==99] <- NA

s2$MC2 <- NA
s2$MC2[s2$Q18_1==1 & s2$Q18_2==4] <- 1
s2$MC2[(s2$Q18_1!=1 | s2$Q18_2!=4) & (!is.na(s2$Q18_1) & !is.na(s2$Q18_2))] <- 0
s2$MC2 <- ifelse(is.na(s2$MC2), 0, s2$MC2)
table(s2$MC2, exclude = NULL)

s2$Q19_1[s2$Q19_1==99] <- NA
s2$Q19_2[s2$Q19_2==99] <- NA

s2$MC3 <- NA
s2$MC3[s2$Q19_1==1 & s2$Q19_2==4] <- 1
s2$MC3[(s2$Q19_1!=1 | s2$Q19_2!=4) & (!is.na(s2$Q19_1) & !is.na(s2$Q19_2))] <- 0
s2$MC3 <- ifelse(is.na(s2$MC3), 0, s2$MC3)
table(s2$MC3, exclude = NULL)

s2$Manipulation_Check <- s2$MC1 + s2$MC2+ s2$MC3
table(s2$Manipulation_Check, s2$Group, exclude = NULL)

s2$Treatment1 <- NA
s2$Treatment1[s2$Group == "Treatment 1"] <- 1
s2$Treatment1[s2$Group == "Control"] <- 0
table(s2$Treatment1, exclude = NULL)

s2$Treatment2 <- NA
s2$Treatment2[s2$Group == "Treatment 2"] <- 1
s2$Treatment2[s2$Group == "Control"] <- 0
table(s2$Treatment2, exclude = NULL)

s2$Treatment3 <- NA
s2$Treatment3[s2$Group == "Treatment 3"] <- 1
s2$Treatment3[s2$Group == "Control"] <- 0
table(s2$Treatment3, exclude = NULL)

bartlett.test(s2$Manipulation_Check ~ s2$Treatment1)
t.test(s2$Manipulation_Check ~ s2$Treatment1, var.equal=T, alternative="less")

bartlett.test(s2$Manipulation_Check ~ s2$Treatment2)
t.test(s2$Manipulation_Check ~ s2$Treatment2, var.equal=T, alternative="less")

bartlett.test(s2$Manipulation_Check ~ s2$Treatment3)
t.test(s2$Manipulation_Check ~ s2$Treatment3, var.equal=F, alternative="less")

#### Supplementary Appendix, Table 14. Study 2 descriptive statistics. ####

library(table1)

my.render.cont <- function(x) {
  with(stats.apply.rounding(stats.default(x), digits=2), c("",
                                                           "Mean (SD)"=sprintf("%s (&plusmn; %s)", MEAN, SD)))
}
my.render.cat <- function(x) {
  c("", sapply(stats.default(x), function(y) with(y,
                                                  sprintf("%d (%0.0f %%)", FREQ, PCT))))
}

table1(~ Self_Conservatism + priming_fear + sex + age + education + self_left_right | Group, data=s2, 
       render.continuous=my.render.cont, render.categorical=my.render.cat, overall=NULL)

treatment <- aov(Self_Conservatism ~ factor(Group), data = s2)
summary(treatment)
TukeyHSD(treatment)

sex <- aov(sex ~ factor(Group), data = s2)
summary(sex)

age <- aov(age ~ factor(Group), data = s2)
summary(age)

education <- aov(education ~ factor(Group), data = s2)
summary(education)

self_left_right <- aov(self_left_right ~ factor(Group), data = s2)
summary(self_left_right)

#### Standardization ####

s2 <- subset(s2, subset = !is.na(Self_Conservatism))

mean(s2$Cancer_Prime, na.rm = T)
s2$beta_Cancer_Prime <- (s2$Cancer_Prime - mean(s2$Cancer_Prime, na.rm = T))
table(s2$beta_Cancer_Prime, exclude = NULL)

mean(s2$Terror_Prime, na.rm = T)
s2$beta_Terror_Prime <- (s2$Terror_Prime - mean(s2$Terror_Prime, na.rm = T))
table(s2$beta_Terror_Prime, exclude = NULL)

mean(s2$General_Mortality_Prime, na.rm = T)
s2$beta_General_Mortality_Prime <- (s2$General_Mortality_Prime - mean(s2$General_Mortality_Prime, na.rm = T))
table(s2$beta_General_Mortality_Prime, exclude = NULL)

mean(s2$sex, na.rm = T)
s2$beta_sex <- (s2$sex - mean(s2$sex, na.rm = T))
table(s2$beta_sex, exclude = NULL)

mean(s2$age, na.rm = T)
sd(s2$age, na.rm = T)
s2$beta_age <- (s2$age - mean(s2$age, na.rm = T))/(2*sd(s2$age, na.rm = T))
table(s2$beta_age, exclude = NULL)

mean(s2$education, na.rm = T)
sd(s2$education, na.rm = T)
s2$beta_education <- (s2$education - mean(s2$education, na.rm = T))/(2*sd(s2$education, na.rm = T))
table(s2$beta_education, exclude = NULL)

mean(s2$self_left_right, na.rm = T)
sd(s2$self_left_right, na.rm = T)
s2$beta_self_left_right <- (s2$self_left_right - mean(s2$self_left_right, na.rm = T))/(2*sd(s2$self_left_right, na.rm = T))
table(s2$beta_self_left_right, exclude = NULL)

mean(s2$Self_Conservatism, na.rm = T)
sd(s2$Self_Conservatism, na.rm = T)
s2$beta_Self_Conservatism <- (s2$Self_Conservatism - mean(s2$Self_Conservatism, na.rm = T))/(2*sd(s2$Self_Conservatism, na.rm = T))
table(s2$beta_Self_Conservatism, exclude = NULL)

mean(s2$priming_fear, na.rm = T)
sd(s2$priming_fear, na.rm = T)
s2$beta_priming_fear <- (s2$priming_fear - mean(s2$priming_fear, na.rm = T))/(2*sd(s2$priming_fear, na.rm = T))
table(s2$beta_priming_fear, exclude = NULL)

s2$State <- ifelse(s2$State=="canakkale", "Canakkale", s2$State)
s2$State <- ifelse(s2$State=="cankiri", "Cankiri", s2$State)
s2$State <- ifelse(s2$State=="kahramanmaras", "Kahramanmaras", s2$State)
s2$State <- ifelse(s2$State=="Kinkkale", "Kirikkale", s2$State)
s2$State <- ifelse(s2$State=="zonguldak", "Zonguldak", s2$State)

s2$region <- NA
s2$region[s2$State=='Istanbul'] <- "Istanbul"
s2$region[s2$State=='Tekirdag' |
            s2$State=='Edirne' |
            s2$State== 'Kirklareli' |
            s2$State== 'Balikesir' | 
            s2$State== 'Canakkale'] <- "West Marmara"
s2$region[s2$State=='Izmir' |
            s2$State== 'Aydin' | 
            s2$State==  'Denizli' |
            s2$State==  'Mugla' |
            s2$State==  'Manisa' |
            s2$State==  'Afyon' |
            s2$State== 'Kutahya' |
            s2$State== 'Usak'] <- "Aegean"
s2$region[s2$State=='Bursa' |
            s2$State== 'Eskisehir' |
            s2$State== 'Bilecik' |
            s2$State== 'Kocaeli' |
            s2$State== 'Sakarya' |
            s2$State==  'Duzce' |
            s2$State==  'Bolu' |
            s2$State== 'Yalova'] <- "East Marmara"
s2$region[s2$State=='Ankara' |
            s2$State==  'Konya' | s2$State== 'Karaman'] <- "West Anatolia"
s2$region[s2$State=='Antalya' |
            s2$State== 'Isparta' |
            s2$State==  'Adana' |
            s2$State==   'Mersin' |
            s2$State==   'Hatay' |
            s2$State==  'Kahramanmaras' |
            s2$State==   'Osmaniye' | s2$State== 'Burdur'] <- "Mediterranean"
s2$region[s2$State=='Kirikkale' |
            s2$State== 'Aksaray' |
            s2$State==  'Nigde' |
            s2$State==   'Nevsehir' |
            s2$State==   'Kirsehir' |
            s2$State==   'Kayseri' | 
            s2$State==  'Sivas' |
            s2$State==  'Yozgat'] <- "Central Anatolia"
s2$region[s2$State=='Zonguldak' |
            s2$State== 'Karabuk' |
            s2$State==  'Bartin' |
            s2$State==  'Kastamonu' |
            s2$State== 'Cankiri' |
            s2$State==  'Sinop' |
            s2$State==  'Samsun' |
            s2$State==  'Tokat' |
            s2$State==  'Corum' |
            s2$State== 'Amasya'] <- "West Black Sea"
s2$region[s2$State=='Ordu' |
            s2$State=='Trabzon' |
            s2$State=='Giresun' |
            s2$State== 'Rize' |
            s2$State== 'Artvin' |
            s2$State== 'Gumushane'] <- "East Black Sea"
s2$region[s2$State=='Erzurum' |
            s2$State=='Erzincan' |
            s2$State== 'Agri' |
            s2$State== 'Kars' | s2$State== 'Ardahan' | s2$State== 'Bayburt' | s2$State== 'Igdir'] <- "Northeast Anatolia"
s2$region[s2$State=='Malatya' |
            s2$State== 'Elazig' |
            s2$State== 'Bingol' |
            s2$State== 'Tunceli' |
            s2$State== 'Van' |
            s2$State== 'Mus' |
            s2$State== 'Bitlis' |
            s2$State== 'Hakkari'] <- "Central East Anatolia"
s2$region[s2$State=='Gaziantep' |
            s2$State== 'Adiyaman' |
            s2$State== 'Kilis' |
            s2$State== 'Sanliurfa' |
            s2$State== 'Diyarbakir' |
            s2$State== 'Mardin' |
            s2$State== 'Batman' |
            s2$State== 'Sirnak' |
            s2$State== 'Siirt'] <- "Southeast Anatolia"
table(s2$region, exclude = NULL)

s2$region <- ifelse(is.na(s2$region), "Outside Turkey", s2$region)

table(s2$State, s2$region)

#### Supplementary Appendix, Table 17. Study 2 balance testing. ####

library(nnet)

balance_test <-  multinom(factor(Group) ~ beta_sex + beta_age + beta_education + beta_self_left_right + beta_akp + beta_chp + beta_other_right + beta_other_left, data=s2)
summary(balance_test)
lmtest::lrtest(balance_test)
tab_model(balance_test, show.ci = 0.95)

p.adjust(c(0.050, 0.711, 0.177, 0.884, 0.199, 0.269, 0.617, 0.680, 0.683), method = p.adjust.methods, n = length(c(0.050, 0.711, 0.177, 0.884, 0.199, 0.269, 0.617, 0.680, 0.683)))
p.adjust(c(0.038, 0.682, 0.095, 0.527, 0.488, 0.357, 0.384, 0.409, 0.868), method = p.adjust.methods, n = length(c(0.038, 0.682, 0.095, 0.527, 0.488, 0.357, 0.384, 0.409, 0.868)))

stargazer(balance_test, type="html", out="s2_balance_test.htm",
          star.cutoffs = c(0.05, 0.01, 0.001),
          digits=2,
          covariate.labels =c("Sex: Male", "Age", "Education", "Left-Right Ideological Placement", "Partisanship: AKP", "Partisanship: CHP", "Partisanship: Other Parties on the Right", "Partisanship: Other Parties on the Left", "Intercept"),
          dep.var.labels   = c("Cancer MS", "Terrorism MS", "General MS"))


#### Supplementary Appendix, Table 20. Study 2 causal mediation analysis models. ####

#### Model 1 ####

out.fit.m01 <- lm(beta_Self_Conservatism ~ beta_Cancer_Prime + beta_Terror_Prime + beta_General_Mortality_Prime, data=s2)
summary(out.fit.m01)
tab_model(out.fit.m01, show.ci = 0.95)

confint(out.fit.m01, 'beta_Cancer_Prime', level=0.95)
confint(out.fit.m01, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m01, 'beta_General_Mortality_Prime', level=0.95)

C <- diag(4)
model.mc <- glht(out.fit.m01, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

med.fit.m1 <- lm(beta_priming_fear ~ beta_Cancer_Prime + beta_Terror_Prime + beta_General_Mortality_Prime, data=s2)
summary(med.fit.m1)
tab_model(med.fit.m1, show.ci = 0.95)

C <- diag(4)
model.mc <- glht(med.fit.m1, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1 <- lm(beta_Self_Conservatism ~ beta_Cancer_Prime + beta_Terror_Prime + beta_General_Mortality_Prime + beta_priming_fear, data=s2)
summary(out.fit.m1)
tab_model(out.fit.m1, show.ci = 0.95)

confint(out.fit.m1, 'beta_Cancer_Prime', level=0.95)
confint(out.fit.m1, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m1, 'beta_General_Mortality_Prime', level=0.95)

C <- diag(5)
model.mc <- glht(out.fit.m1, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

library(mediation)

set.seed(123)
med.out.m1.cancer <- mediate(med.fit.m1, out.fit.m1, treat = "beta_Cancer_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m1.cancer)
med.out.m1.cancer$d.avg.ci

psych::describe(s2$beta_Cancer_Prime)
psych::describe(s2$beta_priming_fear)
psych::describe(s2$beta_Self_Conservatism)

set.seed(123)
med.out.m1.terror <- mediate(med.fit.m1, out.fit.m1, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m1.terror)
med.out.m1.terror$d.avg.ci

psych::describe(s2$beta_Terror_Prime)
psych::describe(s2$beta_priming_fear)
psych::describe(s2$beta_Self_Conservatism)

set.seed(123)
med.out.m1.general <- mediate(med.fit.m1, out.fit.m1, treat = "beta_General_Mortality_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m1.general)
med.out.m1.general$d.avg.ci

psych::describe(s2$beta_General_Mortality_Prime)
psych::describe(s2$beta_priming_fear)
psych::describe(s2$beta_Self_Conservatism)

#### Model 2 ####

out.fit.m02 <- lm(beta_Self_Conservatism ~ beta_Cancer_Prime + beta_Terror_Prime + beta_General_Mortality_Prime + beta_sex + beta_age + beta_education + beta_self_left_right, data=s2)
summary(out.fit.m02)
tab_model(out.fit.m02, show.ci = 0.95)

confint(out.fit.m02, 'beta_Cancer_Prime', level=0.95)
confint(out.fit.m02, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m02, 'beta_General_Mortality_Prime', level=0.95)

C <- diag(8)
model.mc <- glht(out.fit.m02, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

med.fit.m2 <- lm(beta_priming_fear ~ beta_Cancer_Prime + beta_Terror_Prime + beta_General_Mortality_Prime + beta_sex + beta_age + beta_education + beta_self_left_right, data=s2)
summary(med.fit.m2)
tab_model(med.fit.m2, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m2, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m2 <- lm(beta_Self_Conservatism ~ beta_Cancer_Prime + beta_Terror_Prime + beta_General_Mortality_Prime + beta_priming_fear + beta_sex + beta_age + beta_education + beta_self_left_right, data=s2)
summary(out.fit.m2)
tab_model(out.fit.m2, show.ci = 0.95)

confint(out.fit.m2, 'beta_Cancer_Prime', level=0.95)
confint(out.fit.m2, 'beta_Terror_Prime', level=0.95)
confint(out.fit.m2, 'beta_General_Mortality_Prime', level=0.95)

C <- diag(9)
model.mc <- glht(out.fit.m2, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

library(mediation)

set.seed(123)
med.out.m2.cancer <- mediate(med.fit.m2, out.fit.m2, treat = "beta_Cancer_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m2.cancer)

set.seed(123)
med.out.m2.terror <- mediate(med.fit.m2, out.fit.m2, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m2.terror)

set.seed(123)
med.out.m2.general <- mediate(med.fit.m2, out.fit.m2, treat = "beta_General_Mortality_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.out.m2.general)

library(stargazer)

stargazer(out.fit.m01, med.fit.m1, out.fit.m1, out.fit.m02, med.fit.m2, out.fit.m2, 
          type = "html", title=" ", digits=2, out="s2_model_outputs.htm",
          star.cutoffs = c(0.05, 0.01, 0.001),
          model.numbers = F,
          column.labels = c("Model 0.1.1", "Model 1.1", "Model 1.2", "Model 0.2.1", "Model 2.1", "Model 2.2"),
          covariate.labels = c("Cancer Mortality Salience", 
            "Terrorism Mortality Salience", 
            "General Mortality Salience", 
                               "Mortality Fear",
                               "Sex: Male",
                               "Age", 
                               "Education", 
                               "Left-Right Ideological Placement"),
          add.lines = list(c("Estimation", "OLS", "OLS", "OLS", "OLS", "OLS", "OLS")))

#### Supplementary Appendix, Table 21. Study 2 causal moderated mediation analysis models. ####
med.fit.m1.sex <- lm(beta_priming_fear ~ beta_sex*beta_Cancer_Prime + beta_sex*beta_Terror_Prime + beta_sex*beta_General_Mortality_Prime, data=s2)
summary(med.fit.m1.sex)
tab_model(med.fit.m1.sex, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m1.sex, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1.sex <- lm(beta_Self_Conservatism ~ beta_sex*beta_Cancer_Prime + beta_sex*beta_Terror_Prime + beta_sex*beta_General_Mortality_Prime + beta_sex*beta_priming_fear, data=s2)
summary(out.fit.m1.sex)
tab_model(out.fit.m1.sex, show.ci = 0.95)

C <- diag(10)
model.mc <- glht(out.fit.m1.sex, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
round(BH.OUT$test$pvalues, 3)

library(mediation)

table(s2$beta_sex)

set.seed(123)
med.fit.m1.cancer.sex <- mediate(med.fit.m1.sex, out.fit.m1.sex, treat = "beta_Cancer_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.cancer.sex)
test.modmed(med.fit.m1.cancer.sex, covariates.1 = list(beta_sex = -0.336326530612245), covariates.2 = list(beta_sex = 0.663673469387755), sims = 10000)

set.seed(123)
med.fit.m1.terror.sex <- mediate(med.fit.m1.sex, out.fit.m1.sex, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.terror.sex)
test.modmed(med.fit.m1.terror.sex, covariates.1 = list(beta_sex = -0.336326530612245), covariates.2 = list(beta_sex = 0.663673469387755), sims = 10000)

set.seed(123)
med.fit.m1.general.sex <- mediate(med.fit.m1.sex, out.fit.m1.sex, treat = "beta_General_Mortality_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.general.sex)
test.modmed(med.fit.m1.general.sex, covariates.1 = list(beta_sex = -0.336326530612245), covariates.2 = list(beta_sex = 0.663673469387755), sims = 10000)

med.fit.m2.sex <- lm(beta_priming_fear ~ beta_sex*beta_Cancer_Prime + beta_sex*beta_Terror_Prime + beta_sex*beta_General_Mortality_Prime + beta_age + beta_education + beta_self_left_right, data=s2)
summary(med.fit.m2.sex)
tab_model(med.fit.m2.sex, show.ci = 0.95)

C <- diag(11)
model.mc <- glht(med.fit.m2.sex, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m2.sex <- lm(beta_Self_Conservatism ~ beta_sex*beta_Cancer_Prime + beta_sex*beta_Terror_Prime + beta_sex*beta_General_Mortality_Prime + beta_sex*beta_priming_fear + beta_age + beta_education + beta_self_left_right, data=s2)
summary(out.fit.m2.sex)
tab_model(out.fit.m2.sex, show.ci = 0.95)

C <- diag(13)
model.mc <- glht(out.fit.m2.sex, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

set.seed(123)
med.fit.m2.terror.sex <- mediate(med.fit.m2.sex, out.fit.m2.sex, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m2.terror.sex)
test.modmed(med.fit.m2.terror.sex, covariates.1 = list(beta_sex = -0.336326530612245), covariates.2 = list(beta_sex = 0.663673469387755), sims = 10000)

stargazer(med.fit.m1.sex, out.fit.m1.sex, med.fit.m2.sex, out.fit.m2.sex, 
          type = "html", title=" ", digits=2, out="s2_moderated_mediation_outputs.htm",
          star.cutoffs = c(0.05, 0.01, 0.001),
          model.numbers = F,
          column.labels = c("Model 1.1", "Model 1.2", "Model 2.1", "Model 2.2"),
          covariate.labels = c("Sex: Male",
            "Cancer Mortality Salience", 
                               "Terrorism Mortality Salience", 
                               "General Mortality Salience", 
                               "Mortality Fear",
                               "Age", 
                               "Education", 
                               "Left-Right Ideological Placement",
            "Sex: Male x Cancer Mortality Salience", 
            "Sex: Male x Terrorism Mortality Salience", 
            "Sex: Male x General Mortality Salience", 
            "Sex: Male x Mortality Fear"),
          add.lines = list(c("Estimation", "OLS", "OLS", "OLS", "OLS")))

#### Unreported Moderated Mediation by Age ####
med.fit.m1.age <- lm(beta_priming_fear ~ beta_age*beta_Cancer_Prime + beta_age*beta_Terror_Prime + beta_age*beta_General_Mortality_Prime, data=s2)
summary(med.fit.m1.age)
tab_model(med.fit.m1.age, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m1.age, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1.age <- lm(beta_Self_Conservatism ~ beta_age*beta_Cancer_Prime + beta_age*beta_Terror_Prime + beta_age*beta_General_Mortality_Prime + beta_age*beta_priming_fear, data=s2)
summary(out.fit.m1.age)
tab_model(out.fit.m1.age, show.ci = 0.95)

C <- diag(10)
model.mc <- glht(out.fit.m1.age, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

table(s2$beta_age, exclude = NULL)

library(mediation)

set.seed(123)
med.fit.m1.cancer.age <- mediate(med.fit.m1.age, out.fit.m1.age, treat = "beta_Cancer_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.cancer.age)
test.modmed(med.fit.m1.cancer.age, covariates.1 = list(beta_age = -0.334139870797125), covariates.2 = list(beta_age = 3.86353377446537), sims = 10000)

set.seed(123)
med.fit.m1.terror.age <- mediate(med.fit.m1.age, out.fit.m1.age, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.terror.age)
test.modmed(med.fit.m1.terror.age, covariates.1 = list(beta_age = -0.334139870797125), covariates.2 = list(beta_age = 3.86353377446537), sims = 10000)

set.seed(123)
med.fit.m1.general.age <- mediate(med.fit.m1.age, out.fit.m1.age, treat = "beta_General_Mortality_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.general.age)
test.modmed(med.fit.m1.general.age, covariates.1 = list(beta_age = -0.334139870797125), covariates.2 = list(beta_age = 3.86353377446537), sims = 10000)

#### Unreported Moderated Mediation by Left-Right Ideological Placement ####
med.fit.m1.ideology <- lm(beta_priming_fear ~ beta_self_left_right*beta_Cancer_Prime + beta_self_left_right*beta_Terror_Prime + beta_self_left_right*beta_General_Mortality_Prime, data=s2)
summary(med.fit.m1.ideology)
tab_model(med.fit.m1.ideology, show.ci = 0.95)

C <- diag(8)
model.mc <- glht(med.fit.m1.ideology, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

out.fit.m1.ideology <- lm(beta_Self_Conservatism ~ beta_self_left_right*beta_Cancer_Prime + beta_self_left_right*beta_Terror_Prime + beta_self_left_right*beta_General_Mortality_Prime + beta_self_left_right*beta_priming_fear, data=s2)
summary(out.fit.m1.ideology)
tab_model(out.fit.m1.ideology, show.ci = 0.95)

C <- diag(10)
model.mc <- glht(out.fit.m1.ideology, linfct = C)
summary(model.mc)
BH.OUT <- summary(model.mc, test = adjusted(type = "BH"))
BH.OUT

table(s2$beta_self_left_right, exclude = NULL)

library(mediation)

set.seed(123)
med.fit.m1.cancer.ideology <- mediate(med.fit.m1.ideology, out.fit.m1.ideology, treat = "beta_Cancer_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.cancer.ideology)
test.modmed(med.fit.m1.cancer.ideology, covariates.1 = list(beta_self_left_right = -0.770510489116555), covariates.2 = list(beta_self_left_right = 0.78085291178926), sims = 10000)

set.seed(123)
med.fit.m1.terror.ideology <- mediate(med.fit.m1.ideology, out.fit.m1.ideology, treat = "beta_Terror_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.terror.ideology)
test.modmed(med.fit.m1.terror.ideology, covariates.1 = list(beta_self_left_right = -0.770510489116555), covariates.2 = list(beta_self_left_right = 0.78085291178926), sims = 10000)

set.seed(123)
med.fit.m1.general.ideology <- mediate(med.fit.m1.ideology, out.fit.m1.ideology, treat = "beta_General_Mortality_Prime", mediator = "beta_priming_fear", sims = 10000, boot = T, conf.level = 0.95)
summary(med.fit.m1.general.ideology)
test.modmed(med.fit.m1.general.ideology, covariates.1 = list(beta_self_left_right = -0.770510489116555), covariates.2 = list(beta_self_left_right = 0.78085291178926), sims = 10000)

#### Unreported Serial Mediation Analyses (provided for the interested readers) ####
##### Note to the Readers #####

# Serial mediation analyses, testing the alternative causal chain running from treatments to mortality fear to personality traits to conservatism, 
# displayed that the theorized mechanism was not confounded by alternative paths. 
# There was one unpredicted but the significant and nearly replicable alternate mechanism. 
# The participants under General MS got afraid of mortality in general, which decreased perception of seeing the self as conscientious and 
# then led to lesser conservative self-placement (0.01 standard deviations decrease in conservatism with 95% CI [-0.02, -0.001] and 77% statistical power).
# I argue that this path should not be interpreted as abrupt personality changes in the face of threats but rather as adjusting one’s self-perception to alleviate mortality fear (at the moment).

s2$Q20_1[s2$Q20_1==99] <- NA
s2$Q20_2[s2$Q20_2==99] <- NA
s2$Q20_3[s2$Q20_3==99] <- NA
s2$Q20_4[s2$Q20_4==99] <- NA
s2$Q21_1[s2$Q21_1==99] <- NA
s2$Q21_2[s2$Q21_2==99] <- NA
s2$Q21_3[s2$Q21_3==99] <- NA
s2$Q21_4[s2$Q21_4==99] <- NA
s2$Q22_1[s2$Q22_1==99] <- NA
s2$Q22_2[s2$Q22_2==99] <- NA
s2$Q22_3[s2$Q22_3==99] <- NA
s2$Q22_4[s2$Q22_4==99] <- NA
s2$Q23_1[s2$Q23_1==99] <- NA
s2$Q23_2[s2$Q23_2==99] <- NA
s2$Q23_3[s2$Q23_3==99] <- NA
s2$Q23_4[s2$Q23_4==99] <- NA
s2$Q24_1[s2$Q24_1==99] <- NA
s2$Q24_2[s2$Q24_2==99] <- NA
s2$Q24_3[s2$Q24_3==99] <- NA
s2$Q24_4[s2$Q24_4==99] <- NA
s2$Q25_1[s2$Q25_1==99] <- NA
s2$Q25_2[s2$Q25_2==99] <- NA
s2$Q25_3[s2$Q25_3==99] <- NA
s2$Q25_4[s2$Q25_4==99] <- NA
s2$Q26_1[s2$Q26_1==99] <- NA
s2$Q26_2[s2$Q26_2==99] <- NA
s2$Q26_3[s2$Q26_3==99] <- NA
s2$Q26_4[s2$Q26_4==99] <- NA
s2$Q27_1[s2$Q27_1==99] <- NA
s2$Q27_2[s2$Q27_2==99] <- NA
s2$Q27_3[s2$Q27_3==99] <- NA
s2$Q27_4[s2$Q27_4==99] <- NA
s2$Q28_1[s2$Q28_1==99] <- NA
s2$Q28_2[s2$Q28_2==99] <- NA
s2$Q28_3[s2$Q28_3==99] <- NA
s2$Q28_4[s2$Q28_4==99] <- NA
s2$Q29_1[s2$Q29_1==99] <- NA
s2$Q29_2[s2$Q29_2==99] <- NA
s2$Q29_3[s2$Q29_3==99] <- NA
s2$Q29_4[s2$Q29_4==99] <- NA
s2$Q30_1[s2$Q30_1==99] <- NA
s2$Q30_2[s2$Q30_2==99] <- NA
s2$Q30_3[s2$Q30_3==99] <- NA
s2$Q30_4[s2$Q30_4==99] <- NA

names(s2)[names(s2)=="Q20_1"] <- "extraversion_1"
s2$Q21_2R <- 12 - s2$Q21_2
names(s2)[names(s2)=="Q21_2R"] <- "extraversion_2"
names(s2)[names(s2)=="Q22_3"] <- "extraversion_3"
names(s2)[names(s2)=="Q23_4"] <- "extraversion_4"
s2$Q25_1R <- 12 - s2$Q25_1
names(s2)[names(s2)=="Q25_1R"] <- "extraversion_5"
names(s2)[names(s2)=="Q26_2"] <- "extraversion_6"
s2$Q27_3R <- 12 - s2$Q27_3
names(s2)[names(s2)=="Q27_3R"] <- "extraversion_7"
names(s2)[names(s2)=="Q28_4"] <- "extraversion_8"

psych::alpha(s2[s2$Group=="Control", c("extraversion_1", "extraversion_2", "extraversion_3", "extraversion_4",
                   "extraversion_5", "extraversion_6", "extraversion_7", "extraversion_8")])

s2$extraversion <- s2$extraversion_1 + 
  s2$extraversion_2 + 
  s2$extraversion_3 + 
  s2$extraversion_4 + 
  s2$extraversion_5 + 
  s2$extraversion_6 + 
  s2$extraversion_7 + 
  s2$extraversion_8

psych::describeBy(s2$extraversion, s2$Group)

mean(s2$extraversion, na.rm = T)
sd(s2$extraversion, na.rm = T)
s2$beta_extraversion <- (s2$extraversion - mean(s2$extraversion, na.rm = T))/(2*sd(s2$extraversion, na.rm = T))
table(s2$beta_extraversion, exclude = NULL)

s2$Q20_2R <- 12 - s2$Q20_2
names(s2)[names(s2)=="Q20_2R"] <- "agreeableness_1"
names(s2)[names(s2)=="Q21_3"] <- "agreeableness_2"
s2$Q22_4R <- 12 - s2$Q22_4
names(s2)[names(s2)=="Q22_4R"] <- "agreeableness_3"
names(s2)[names(s2)=="Q24_1"] <- "agreeableness_4"
names(s2)[names(s2)=="Q25_2"] <- "agreeableness_5"
s2$Q26_3R <- 12 - s2$Q26_3
names(s2)[names(s2)=="Q26_3R"] <- "agreeableness_6"
names(s2)[names(s2)=="Q27_4"] <- "agreeableness_7"
s2$Q29_1R <- 12 - s2$Q29_1
names(s2)[names(s2)=="Q29_1R"] <- "agreeableness_8"
names(s2)[names(s2)=="Q30_2"] <- "agreeableness_9"

psych::alpha(s2[s2$Group=="Control", c("agreeableness_1", "agreeableness_2", "agreeableness_3", "agreeableness_4",
                   "agreeableness_5", "agreeableness_6", "agreeableness_7", "agreeableness_8", "agreeableness_9")])

s2$agreeableness <- s2$agreeableness_1 + 
  s2$agreeableness_2 + 
  s2$agreeableness_3 + 
  s2$agreeableness_4 + 
  s2$agreeableness_5 + 
  s2$agreeableness_6 + 
  s2$agreeableness_7 + 
  s2$agreeableness_8 +
  s2$agreeableness_9

psych::describeBy(s2$agreeableness, s2$Group)

mean(s2$agreeableness, na.rm = T)
sd(s2$agreeableness, na.rm = T)
s2$beta_agreeableness <- (s2$agreeableness - mean(s2$agreeableness, na.rm = T))/(2*sd(s2$agreeableness, na.rm = T))
table(s2$beta_agreeableness, exclude = NULL)

names(s2)[names(s2)=="Q20_3"] <- "conscientiousness_1"
s2$Q21_4R <- 12 - s2$Q21_4
names(s2)[names(s2)=="Q21_4R"] <- "conscientiousness_2"
names(s2)[names(s2)=="Q23_1"] <- "conscientiousness_3"
s2$Q24_2R <- 12 - s2$Q24_2
names(s2)[names(s2)=="Q24_2R"] <- "conscientiousness_4"
s2$Q25_3R <- 12 - s2$Q25_3
names(s2)[names(s2)=="Q25_3R"] <- "conscientiousness_5"
names(s2)[names(s2)=="Q26_4"] <- "conscientiousness_6"
names(s2)[names(s2)=="Q28_1"] <- "conscientiousness_7"
names(s2)[names(s2)=="Q29_2"] <- "conscientiousness_8"
s2$Q30_3R <- 12 - s2$Q30_3
names(s2)[names(s2)=="Q30_3R"] <- "conscientiousness_9"

psych::alpha(s2[s2$Group=="Control", c("conscientiousness_1", "conscientiousness_2", "conscientiousness_3", "conscientiousness_4",
                   "conscientiousness_5", "conscientiousness_6", "conscientiousness_7", "conscientiousness_8", "conscientiousness_9")])

s2$conscientiousness <- s2$conscientiousness_1 + 
  s2$conscientiousness_2 + 
  s2$conscientiousness_3 + 
  s2$conscientiousness_4 + 
  s2$conscientiousness_5 + 
  s2$conscientiousness_6 + 
  s2$conscientiousness_7 + 
  s2$conscientiousness_8 +
  s2$conscientiousness_9

psych::describeBy(s2$conscientiousness, s2$Group)

mean(s2$conscientiousness, na.rm = T)
sd(s2$conscientiousness, na.rm = T)
s2$beta_conscientiousness <- (s2$conscientiousness - mean(s2$conscientiousness, na.rm = T))/(2*sd(s2$conscientiousness, na.rm = T))
table(s2$beta_conscientiousness, exclude = NULL)

names(s2)[names(s2)=="Q20_4"] <- "neuroticism_1"
s2$Q22_1R <- 12 - s2$Q22_1
names(s2)[names(s2)=="Q22_1R"] <- "neuroticism_2"
names(s2)[names(s2)=="Q23_2"] <- "neuroticism_3"
names(s2)[names(s2)=="Q24_3"] <- "neuroticism_4"
s2$Q25_4R <- 12 - s2$Q25_4
names(s2)[names(s2)=="Q25_4R"] <- "neuroticism_5"
names(s2)[names(s2)=="Q27_1"] <- "neuroticism_6"
s2$Q28_2R <- 12 - s2$Q28_2
names(s2)[names(s2)=="Q28_2R"] <- "neuroticism_7"
names(s2)[names(s2)=="Q29_3"] <- "neuroticism_8"

psych::alpha(s2[s2$Group=="Control", c("neuroticism_1", "neuroticism_2", "neuroticism_3", "neuroticism_4",
                   "neuroticism_5", "neuroticism_6", "neuroticism_7", "neuroticism_8")])

s2$neuroticism <- s2$neuroticism_1 + 
  s2$neuroticism_2 + 
  s2$neuroticism_3 + 
  s2$neuroticism_4 + 
  s2$neuroticism_5 + 
  s2$neuroticism_6 + 
  s2$neuroticism_7 + 
  s2$neuroticism_8

psych::describeBy(s2$neuroticism, s2$Group)

mean(s2$neuroticism, na.rm = T)
sd(s2$neuroticism, na.rm = T)
s2$beta_neuroticism <- (s2$neuroticism - mean(s2$neuroticism, na.rm = T))/(2*sd(s2$neuroticism, na.rm = T))
table(s2$beta_neuroticism, exclude = NULL)

names(s2)[names(s2)=="Q21_1"] <- "openness_1"
names(s2)[names(s2)=="Q22_2"] <- "openness_2"
names(s2)[names(s2)=="Q23_3"] <- "openness_3"
names(s2)[names(s2)=="Q24_4"] <- "openness_4"
names(s2)[names(s2)=="Q26_1"] <- "openness_5"
names(s2)[names(s2)=="Q27_2"] <- "openness_6"
s2$Q28_3R <- 12 - s2$Q28_3
names(s2)[names(s2)=="Q28_3R"] <- "openness_7"
names(s2)[names(s2)=="Q29_4"] <- "openness_8"
s2$Q30_1R <- 12 - s2$Q30_1
names(s2)[names(s2)=="Q30_1R"] <- "openness_9"
names(s2)[names(s2)=="Q30_4"] <- "openness_10"

psych::alpha(s2[s2$Group=="Control", c("openness_1", "openness_2", "openness_3", "openness_4",
                   "openness_5", "openness_6", "openness_7", "openness_8", "openness_9", "openness_10")])

s2$openness <- s2$openness_1 + 
  s2$openness_2 + 
  s2$openness_3 + 
  s2$openness_4 + 
  s2$openness_5 + 
  s2$openness_6 + 
  s2$openness_7 + 
  s2$openness_8 +
  s2$openness_9 +
  s2$openness_10

psych::describeBy(s2$openness, s2$Group)

mean(s2$openness, na.rm = T)
sd(s2$openness, na.rm = T)
s2$beta_openness <- (s2$openness - mean(s2$openness, na.rm = T))/(2*sd(s2$openness, na.rm = T))
table(s2$beta_openness, exclude = NULL)

#### Unreported Study 2 descriptive statistics of the personality traits. ####

table1(~ openness + conscientiousness + extraversion + agreeableness + neuroticism | Group, data=s2, 
       render.continuous=my.render.cont, render.categorical=my.render.cat, overall=NULL)

library(lavaan)

base.model.1.acme =
  "
  beta_priming_fear ~ c_f*beta_Cancer_Prime + t_f*beta_Terror_Prime + d_f*beta_General_Mortality_Prime
  
  beta_openness ~ c_b5op*beta_Cancer_Prime + t_b5op*beta_Terror_Prime + d_b5op*beta_General_Mortality_Prime + f_b5op*beta_priming_fear + beta_sex + beta_age + beta_education + beta_self_left_right
  beta_conscientiousness ~ c_b5co*beta_Cancer_Prime + t_b5co*beta_Terror_Prime + d_b5co*beta_General_Mortality_Prime + f_b5co*beta_priming_fear + beta_sex + beta_age + beta_education + beta_self_left_right
  beta_extraversion ~ c_b5ex*beta_Cancer_Prime + t_b5ex*beta_Terror_Prime + d_b5ex*beta_General_Mortality_Prime + f_b5ex*beta_priming_fear + beta_sex + beta_age + beta_education + beta_self_left_right
  beta_agreeableness ~ c_b5ag*beta_Cancer_Prime + t_b5ag*beta_Terror_Prime + d_b5ag*beta_General_Mortality_Prime + f_b5ag*beta_priming_fear + beta_sex + beta_age + beta_education + beta_self_left_right
  beta_neuroticism ~ c_b5ne*beta_Cancer_Prime + t_b5ne*beta_Terror_Prime + d_b5ne*beta_General_Mortality_Prime + f_b5ne*beta_priming_fear + beta_sex + beta_age + beta_education + beta_self_left_right
  
  beta_Self_Conservatism ~ c_con*beta_Cancer_Prime + t_con*beta_Terror_Prime + d_con*beta_General_Mortality_Prime + f_con*beta_priming_fear +
  b5op_con*beta_openness + b5co_con*beta_conscientiousness + b5ex_con*beta_extraversion + b5ag_con*beta_agreeableness + b5ne_con*beta_neuroticism +
  beta_sex + beta_age + beta_education + beta_self_left_right
  
  acme_cancer_fear := c_f*f_con
  acme_terrorism_fear := t_f*f_con
  acme_General_Mortality_fear := d_f*f_con
  
  acme_cancer_fear_b5op := c_f*f_b5op*b5op_con
  acme_cancer_fear_b5co := c_f*f_b5co*b5co_con
  acme_cancer_fear_b5ex := c_f*f_b5ex*b5ex_con
  acme_cancer_fear_b5ag := c_f*f_b5ag*b5ag_con
  acme_cancer_fear_b5ne := c_f*f_b5ne*b5ne_con
  
  acme_terrorism_fear_b5op := t_f*f_b5op*b5op_con
  acme_terrorism_fear_b5co := t_f*f_b5co*b5co_con
  acme_terrorism_fear_b5ex := t_f*f_b5ex*b5ex_con
  acme_terrorism_fear_b5ag := t_f*f_b5ag*b5ag_con
  acme_terrorism_fear_b5ne := t_f*f_b5ne*b5ne_con
  
  acme_General_Mortality_fear_b5op := d_f*f_b5op*b5op_con
  acme_General_Mortality_fear_b5co := d_f*f_b5co*b5co_con
  acme_General_Mortality_fear_b5ex := d_f*f_b5ex*b5ex_con
  acme_General_Mortality_fear_b5ag := d_f*f_b5ag*b5ag_con
  acme_General_Mortality_fear_b5ne := d_f*f_b5ne*b5ne_con
  
  acme_cancer_b5co := c_b5ne*b5co_con
  acme_terrorism_b5co := t_b5ne*b5co_con
  acme_General_Mortality_b5co := d_b5co*b5ne_con
  
  acme_cancer_b5ne := c_b5ne*b5ne_con
  acme_terrorism_b5ne := t_b5ne*b5ne_con
  acme_General_Mortality_b5ne := d_b5ne*b5ne_con
"
fit.base.model.1.acme <-  sem(base.model.1.acme, data=s2, se = "bootstrap", bootstrap = 10000)
summary(fit.base.model.1.acme, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE, nd = 2)
parameterEstimates(fit.base.model.1.acme, rsquare=T, output = "text", level = 0.95)

psych::describe(s2$beta_Cancer_Prime)
psych::describe(s2$beta_Terror_Prime)
psych::describe(s2$beta_General_Mortality_Prime)

psych::describe(s2$beta_priming_fear)

psych::describe(s2$beta_conscientiousness)
psych::describe(s2$beta_neuroticism)

psych::describe(s2$beta_Self_Conservatism)

#### Supplementary Appendix, Study 2 Sample Representativeness Tests, Tables 5 to 8 ####

table(s2$age, exclude = NULL)
s2$age_cat <- NA
s2$age_cat[s2$age>=18 & s2$age<=24] <- "18-24"
s2$age_cat[s2$age>=25 & s2$age<=34] <- "25-34"
s2$age_cat[s2$age>=35 & s2$age<=49] <- "35-49"
s2$age_cat[s2$age>=50 & s2$age<=64] <- "50-64"
s2$age_cat[s2$age>=65] <- "65 or Over"
table(s2$age_cat, exclude = NULL)
round(prop.table(table(s2$age_cat, exclude = NULL))*100, 1)

s2_age <- c(75.1, 16.6, 5.8, 2.1, 0.4)
tr_2018_age <- c(15.3, 21.3, 30.1, 21.1, 12.3)

age_s2_tr_2018 <- cbind(s2_age, tr_2018_age)
row.names(age_s2_tr_2018) <- c("18-24","25-34", "35-49", "50-64","65+")

chisq.test(age_s2_tr_2018)
round(chisq.test(age_s2_tr_2018 )$stdres, 2)
round(rcompanion::cramerV(age_s2_tr_2018), 2)

table(s2$sex, exclude = NULL)
s2$sex_cat <- NA
s2$sex_cat[s2$sex==1] <- "Male"
s2$sex_cat[s2$sex==0] <- "Female"
table(s2$sex_cat, exclude = NULL)
round(prop.table(table(s2$sex_cat, exclude = NULL))*100, 1)

s2_gender <- c(33.1, #male
               65.2) #female

tr_2018_gender <- c(49.2,
                    50.8)

gender_s2_tr_2018 <- cbind(s2_gender, tr_2018_gender)
row.names(gender_s2_tr_2018) <- c("Male","Female")

chisq.test(gender_s2_tr_2018)
round(rcompanion::cramerV(gender_s2_tr_2018), 2)

table(s2$education, exclude = NULL)
s2$education_cat <- NA
s2$education_cat[s2$education==2 | s2$education==3] <- "Below High School Graduate"
s2$education_cat[s2$education==4 | s2$education==5] <- "High School Graduate"
s2$education_cat[s2$education==6 | s2$education==7 | s2$education==8] <- "College Graduate or Above"
s2$education_cat[is.na(s2$education)] <- "Unknown"
table(s2$education_cat, exclude = NULL)
s2$education_cat <- factor(s2$education_cat, levels = c("Below High School Graduate", "High School Graduate",
                                                        "College Graduate or Above", "Unknown"))
round(prop.table(table(s2$education_cat, exclude = NULL))*100, 1)

s2_education <- c(2.5, 69.7, 26.6, 1.2)
tr_2018_education <- c(47.9, 24.2, 18, 0.7)

education_s2_tr_2018 <- cbind(s2_education, tr_2018_education)
row.names(education_s2_tr_2018) <- c("Below High School Graduate",
                                     "High School Graduate",
                                     "College Graduate or Above",
                                     "Unknown")

chisq.test(education_s2_tr_2018)
round(chisq.test(education_s2_tr_2018 )$stdres, 2)
round(rcompanion::cramerV(education_s2_tr_2018), 2)

round(prop.table(table(s2$State, exclude = NULL))*100, 1)

s2_geo <- c(0.1,
            0.6,
            0.9,
            0.2,
            0.2,
            0.9,
            0.4,
            9.2,
            1,
            0.2,
            0.1,
            0,
            0.2,
            0.1,
            0.1,
            0.2,
            0.7,
            0.4,
            0,
            0.1,
            0,
            0.3,
            0.2,
            0.1,
            0,
            0,
            0.5,
            0.1,
            0.3,
            0.2,
            1.7,
            0,
            0.7,
            0.6,
            0.2,
            3.3,
            0.2,
            0.6,
            44.6,
            0.4,
            0.3,
            0.2,
            0,
            1.9,
            0.6,
            2.7,
            3.8,
            0.2,
            0.1,
            0,
            0,
            0,
            0.2,
            0.2,
            0.1,
            0.4,
            1.1,
            0.1,
            0,
            0.5,
            0.2,
            0.1,
            0.1,
            2,
            0,
            10.8,
            0,
            0.2,
            0.8,
            0.1,
            0.6,
            0,
            0.2,
            0.2,
            0.3,
            0.2,
            0.6,
            0.2,
            0.1,
            1,
            0.2)

tr_2018_geo <- c(1.9,
                 1.4,
                 0.9,
                 1.3,
                 0.7,
                 0.5,
                 1.3,
                 5.7,
                 1.7,
                 0.8,
                 0.5,
                 0.4,
                 0.4,
                 0.3,
                 0.4,
                 0.5,
                 1,
                 1.1,
                 0.7,
                 0.3,
                 0.3,
                 0.4,
                 0.4,
                 0.1,
                 1,
                 1,
                 0.6,
                 0.4,
                 0.2,
                 0.2,
                 1.3,
                 0.5,
                 1.1,
                 0.3,
                 0.4,
                 3.8,
                 0.3,
                 2.4,
                 18.7,
                 1.3,
                 1.8,
                 2.3,
                 0.6,
                 0.6,
                 2.7,
                 3,
                 0.3,
                 0.2,
                 0.5,
                 0.1,
                 0.1,
                 0.3,
                 0.9,
                 0.3,
                 1.9,
                 1.8,
                 2.2,
                 0.8,
                 0.5,
                 0.6,
                 0.2,
                 0.7,
                 0.3,
                 2.7,
                 0.3,
                 6.9,
                 1.7,
                 0.8,
                 0.7,
                 0.4,
                 0.3,
                 0.3,
                 0.3,
                 0.2,
                 0.5,
                 0.8,
                 1.7,
                 1.3,
                 0.7,
                 0.6,
                 0.5)

geo_s2_tr_2018 <- cbind(s2_geo, tr_2018_geo)
row.names(geo_s2_tr_2018) <- c(
  "Manisa",
  "Aydin",
  "Afyon",
  "Mugla",
  "Kutahya",
  "Usak",
  "Denizli",
  "Izmir",
  "Kayseri",
  "Sivas",
  "Aksaray",
  "Nevsehir",
  "Nigde",
  "Kirsehir",
  "Kirikkale",
  "Yozgat",
  "Malatya",
  "Van",
  "Elazig",
  "Hakkari",
  "Bingol",
  "Bitlis",
  "Mus",
  "Tunceli",
  "Ordu",
  "Trabzon",
  "Giresun",
  "Rize",
  "Artvin",
  "Gumushane",
  "Sakarya",
  "Duzce",
  "Eskisehir",
  "Bilecik",
  "Bolu",
  "Bursa",
  "Yalova",
  "Kocaeli",
  "Istanbul",
  "Kahramanmaras",
  "Hatay",
  "Mersin",
  "Isparta",
  "Osmaniye",
  "Adana",
  "Antalya",
  "Burdur",
  "Igdir",
  "Agri",
  "Ardahan",
  "Bayburt",
  "Erzincan",
  "Erzurum",
  "Kars",
  "Sanliurfa",
  "Diyarbakir",
  "Gaziantep",
  "Mardin",
  "Sirnak",
  "Batman",
  "Kilis",
  "Adiyaman",
  "Siirt",
  "Konya",
  "Karaman",
  "Ankara",
  "Samsun",
  "Tokat",
  "Corum",
  "Amasya",
  "Bartin",
  "Karabuk",
  "Sinop",
  "Cankiri",
  "Kastamonu",
  "Zonguldak",
  "Balikesir",
  "Tekirdag",
  "Canakkale",
  "Edirne",
  "Kirklareli")

chisq.test(geo_s2_tr_2018)
round(chisq.test(geo_s2_tr_2018)$stdres, 2)
round(rcompanion::cramerV(geo_s2_tr_2018), 2)
